diff --git a/main.F b/main.F
index b735ca9..4fb44a5 100644
--- a/main.F
+++ b/main.F
@@ -1304,6 +1304,8 @@
               WRITE(TIU6,*)'finite differences'
       ELSE IF (DYN%IBRION==10) THEN
            WRITE(TIU6,*)'relaxation of ions and charge simultaneously'
+      ELSE IF (DYN%IBRION==12) THEN
+           WRITE(TIU6,*)'REFTRAJ mode'
       ENDIF
 
       IF (DYN%IBRION/=-1 .AND. T_INFO%LSDYN) THEN
@@ -2984,6 +2986,7 @@
         CALL STEP(DYN%INIT,ISCALE,T_INFO%NIONS,LATT_CUR%A,LATT_CUR%ANORM,DYN%D2C,DYN%SMASS,DYN%POSION,DYN%POSIOC, &
              DYN%POTIM,T_INFO%POMASS,T_INFO%NTYP,T_INFO%ITYP,DYN%TEMP,DYN%VEL,DYN%D2,DYN%D3,DYN%SNOSE, &
              EKIN,EPS,ES,DISMAX,NDEGREES_OF_FREEDOM, IO%IU6)
+
         TEIN = 2*EKIN/BOLKEV/NDEGREES_OF_FREEDOM
 
 ! sum energy of images along chain
@@ -3289,6 +3292,8 @@
 ! if we have very small forces (small trial energy change) we can stop
            IF (IFLAG==1) INFO%LSTOP=INFO%LSTOP .OR. (ABS(E1TEST) < 0.1_q*DYN%EDIFFG)
 !-----------------------------------------------------------------------
+        ELSE IF (DYN%IBRION == 12) THEN
+          CALL REFTRAJ(LATT_CUR, T_INFO, DYN, TOTEN, TSIF, TIFOR, NSTEP, IO%IU0, INFO%LSTOP, WDES%COMM)
         ENDIF
 
 ! restrict volume for constant volume relaxation
@@ -3303,19 +3308,24 @@
            ENDDO
         ENDIF
         CALL LATTIC(LATT_CUR)
+! for IBRION==12, more like  MD IBRION==0
+	IF (DYN%IBRION==12) THEN
+           PRED%INIPRE=0
+        ELSE
 !  reinitialize the prediction algorithm for the wavefunction if needed
-        PRED%INIPRE=3
-        IF ( PRED%IWAVPR >=12 .AND. &
-             &     (ABS(TOTEN-TOTENG)/T_INFO%NIONS>1.0_q .OR. IFLAG==1)) THEN
-           CALL WAVPRE_NOIO(GRIDC,P,PRED,T_INFO,W,WDES,LATT_CUR,IO%LOPEN, &
-                CHTOT,RHOLM,N_MIX_PAW, CSTRF, LMDIM,CQIJ,INFO%LOVERL,NBLK,IO%IU0)
-
-        ELSE IF ( PRED%IWAVPR >=2 .AND. PRED%IWAVPR <10   .AND. &
-             &     (ABS(TOTEN-TOTENG)/T_INFO%NIONS>1.0_q .OR. IFLAG==1)) THEN
-           CALL WAVPRE(GRIDC,P,PRED,T_INFO,W,WDES,LATT_CUR,IO%LOPEN, &
-                CHTOT,RHOLM,N_MIX_PAW, CSTRF, LMDIM,CQIJ,INFO%LOVERL,NBLK,IO%IU0)
-        ENDIF
-
+           PRED%INIPRE=3
+           IF ( PRED%IWAVPR >=12 .AND. &
+                &     (ABS(TOTEN-TOTENG)/T_INFO%NIONS>1.0_q .OR. IFLAG==1)) THEN
+              CALL WAVPRE_NOIO(GRIDC,P,PRED,T_INFO,W,WDES,LATT_CUR,IO%LOPEN, &
+                   CHTOT,RHOLM,N_MIX_PAW, CSTRF, LMDIM,CQIJ,INFO%LOVERL,NBLK,IO%IU0)
+
+           ELSE IF ( PRED%IWAVPR >=2 .AND. PRED%IWAVPR <10   .AND. &
+                &     (ABS(TOTEN-TOTENG)/T_INFO%NIONS>1.0_q .OR. IFLAG==1)) THEN
+              CALL WAVPRE(GRIDC,P,PRED,T_INFO,W,WDES,LATT_CUR,IO%LOPEN, &
+                   CHTOT,RHOLM,N_MIX_PAW, CSTRF, LMDIM,CQIJ,INFO%LOVERL,NBLK,IO%IU0)
+           ENDIF
+        END IF
+        
         ! use forces as stopping criterion if EDIFFG<0
         IF (DYN%EDIFFG<0) INFO%LSTOP=LSTOP2
         io_begin
diff --git a/poscar.F b/poscar.F
index e9f23fc..39d896f 100644
--- a/poscar.F
+++ b/poscar.F
@@ -591,6 +591,95 @@
       RETURN
       END SUBROUTINE
 
+!=======================================================================
+!
+! read positions from REFTRAJCAR file when it's ready
+!
+!=======================================================================
+      SUBROUTINE REFTRAJ(LATT_CUR, T_INFO, DYN, TOTEN, TSIF, TIFOR, NSTEP, IU0, LSTOP, MYCOMM)
+      USE prec
+      USE lattice
+      USE main_mpi
+      IMPLICIT NONE
+
+      TYPE (latt)::       LATT_CUR
+      TYPE (type_info) :: T_INFO
+      TYPE (dynamics)  :: DYN
+      REAL(q) TOTEN, TSIF(3,3), TIFOR(:,:)
+      INTEGER :: NSTEP
+      INTEGER :: IU0
+      LOGICAL :: LSTOP
+      TYPE (communic) :: MYCOMM
+    ! local
+      INTEGER I, NI, IERROR 
+      LOGICAL exists
+
+      LSTOP = .false.
+
+      ! save POSIOC just for fun
+      DYN%POSIOC = DYN%POSION
+
+#if defined(MPI)
+      if (MYCOMM%IONODE == MYCOMM%NODE_ME) then ! on head node, do the I/O
+#endif
+	 open(file="REFTRAJ_OUTPUT", unit=100, status="UNKNOWN")
+	 write(unit=100,fmt=*) T_INFO%NIONS
+	 write(unit=100,fmt='(F25.16)') TOTEN
+	 do I=1, T_INFO%NIONS
+	    write(unit=100,fmt='(3F25.16)') TIFOR(1,I), TIFOR(2,I), TIFOR(3,I)
+	 end do
+	 write(unit=100,fmt='(6F25.16)') TSIF(1,1), TSIF(2,2), TSIF(3,3), TSIF(1,2), TSIF(2,3), TSIF(3,1)
+	 close(unit=100)
+
+	 open(file="REFTRAJ_STEP_DONE", unit=100, status="UNKNOWN")
+	 write(unit=100,fmt=*) NSTEP
+	 close(unit=100)
+
+	 ! wait for REFTRAJ_READY to exist
+	 exists = .false.
+	 do while (.not. exists)
+	    INQUIRE(file="REFTRAJ_READY", EXIST=exists)
+	    call usleep(100000)
+	 end do
+	 call unlink("REFTRAJ_READY"//char(0))
+	 ! read REFTRAJCAR
+	 open (unit=100, FILE="REFTRAJCAR", STATUS="OLD")
+	 read (unit=100, fmt=*) NI
+	 ! check for NI=0, or invalid NI
+	 if (NI /= 0 .and. NI == T_INFO%NIONS) then
+	    read (unit=100, fmt=*) LATT_CUR%A(1,1:3)
+	    read (unit=100, fmt=*) LATT_CUR%A(2,1:3)
+	    read (unit=100, fmt=*) LATT_CUR%A(3,1:3)
+	    do I=1, NI
+	       read (unit=100, fmt=*) DYN%POSION(1:3,I)
+	    end do
+	 endif
+	 close (unit=100)
+#if defined(MPI)
+      else ! other nodes
+	 NI = 0
+	 LATT_CUR%A = 0.0
+	 DYN%POSION = 0.0
+      endif
+      ! coalesce values onto all nodes
+      CALLMPI( M_sum_i(MYCOMM, NI, 1))
+      CALLMPI( M_sum_d(MYCOMM, LATT_CUR%A(1,1), 3*3))
+      CALLMPI( M_sum_d(MYCOMM, DYN%POSION(1,1), T_INFO%NIONS*3))
+#endif
+      if (NI == 0) then
+	 if (IU0 >= 0) write(iu0,'(A)') "REFTRAJ got NIONS = 0, quitting"
+	 LSTOP = .true.
+      else if (NI /= T_INFO%NIONS) then
+	 if (IU0 >= 0) write(iu0,'(A,2I6)') "REFTRAJ got mismatched atom number ", NI, T_INFO%NIONS
+	 stop
+      endif
+
+      ! fix reciprocal lattice and stuff, just in case
+      CALL LATTIC(LATT_CUR)
+
+      END SUBROUTINE REFTRAJ
+
+
 
 !*************************SUBROUTINE OUTPOS_TRAIL  *********************
 ! write trailer for CONTCAR file
-- 
1.7.2.2

